Numpy Basics

In this notebook we introduce the numpy module. It is perhaps the most widely used library for scientific computing in Python. Many other libraries like the pandas library used in data analysis and the SciPy ecosystem are built on top of it.

Perhaps the main feature of the numpy library is the multidimensional array data structure. Contrary to lists and other iterable types we have seen, numpy arrays were specially designed with linear algebra in mind. This means that all the usual operations from linear algebra are already implemented in numpy with performance and practicality in mind. Let us consider our first example to illustrate this.

import numpy as np
myFirstList = [1,2,3,4]
mySecondList = [7,8,9,10]
plusList = myFirstList + mySecondList
myFirstNumpyArray = np.array(myFirstList)
mySecondNumpyArray = np.array(mySecondList)
plusNumpyArray = myFirstNumpyArray + mySecondNumpyArray
print("Using + with lists: {} + {} = {}".format(myFirstList, mySecondList, plusList))
print("Using + with numpy arrays: {} + {} = {}".format(myFirstNumpyArray, mySecondNumpyArray, plusNumpyArray))

In the example above, note the line import numpy as np. This is the usual way people import the numpy library in planet Earth. You should also pay special attention to the line myFirstNumpyArray = np.array(myFirstList). It summarizes the usual way we create numpy arrays: passing a list as argument to the function np.array(), as below.

myArray = np.array(myList)

Exercise 1. Vector addition with numpy.

Run the code above and compare the behavior of the + operation with lists and numpy arrays.

As you have seen, lists have a computer science bias. They store data, and the "+" operator simply means storing more stuff. The numpy array on the other hand impersonates the physicist's favorite tool: the vector. For numpy arrays, "+" is simply vector addition. And indeed, this reasoning goes much further. We can do to numpy arrays pretty much everything we would do to vectors and get the results we intuitively expect.

Exercise 2. Vector multiplication by a scalar.

Multiply the numpy array np.array([1,2,3,4]) by 4.

Exercise 3. Yet more vectorized operations.

Let a = np.array([1,2,3,4]), b = np.array([5,6,7,8]). Compute a*b, 2**a, a**2, a**b

By now you should already have a good intuitive ideia about how numpy arrays work. We can further demonstrate their convenience with more examples from linear algebra. Suppose we want to transpose a matrix, say

Without numpy, we could do something like this.

def transpose(aMatrix):
    """
    Returns the transpose of a matrix

    aMatrix: a list of embedded lists representing a matrix.
    """
    transp = []
    for column in range(len(aMatrix[0])):
        columnVector = []
        for line in range(len(aMatrix)):
            columnVector += [aMatrix[line][column]]
        transp += [columnVector]
    return transp
a = [[1,2], [3,4], [5,6]]
print("This is a: ", a)
print("This is a transposed: ", transpose(a))
transpose(a)

There are many other ways to do this with lists in Python, of course. But we could otherwise use numpy and simply write

a = np.array([[1,2], [3,4], [5,6]])
a.transpose()

Exercise 4. How would you write the transpose function with list comprehensions?

Exercise 5. The dot product.

The dot (or inner) product of two vectors is found virtually everywhere in science. The usual definition in $\mathbb{R}^n$ in terms of orthonormal vector components $u = (u_0,... ,u_{n-1})$ is $u\cdot v = \sum_{i=0}^{n-1} u_iv_i$. This type of operation in which indices are contracted are far more general, though. As another example, the matrix product of two matrices $A = [A_{ij}]$ and $B = [B_{ij}]$ is given by $[A\cdot B]_{ij} = \sum_{k=0}^{n-1}A_{ik}B_{kj}$. Evidently, we assume that the number of columns in $A$ is identical to the number of lines in $B$. In what follows, implement code that computes the dot and matrix products with lists.

Exercise 6. Benchmarking our code against numpy.

The code below compares the performance of our homemade dot_product() against numpy's np.dot() (you can find the documentation on this function here). Run it and then do the same with our homemade matrix_product() against numpy's np.dot(). Use square matrices. Could you have predicted the shapes of the resulting curves in both comparisons?

# We now time the performance of dot_product() against np.dot().
# We sample two arrays from the random uniform distribution in [0,1),
# with sizes ranging from 4 to 10000 in steps of 4.

def benchmark(aFunction, firstArg, secondArg, numberOfTests = 1):
    """
    Returns the best (that is, the least) processing time spent by a function  
    that takes two arguments.

    aFunction: The function whose performance we are measuring. It takes two arguments.
    firstArg: The first argument of aFunction
    secondArg: The second argument of aFunction
    numberOfTests: The number of times we compute 
                   aFunction(firstArg, secondArg)
    """
    timesTaken = []
    for i in range(numberOfTests):
        start = time.time()
        aFunction(firstArg, secondArg)
        end = time.time()
        timesTaken.append(end-start)
    return min(timesTaken)


dotProductTimes = []
npDotTimes = []
for size in range(10, 10000, 10):
    array1 = np.random.uniform(size = size)
    array2 = np.random.uniform(size = size)
    list1 = list(array1)
    list2 = list(array2)
    timeDotProduct =  benchmark(dot_product, list1, list2,numberOfTests=10)
    timeNpDot = benchmark(np.dot, array1, array2, numberOfTests=10)
    dotProductTimes.append(timeDotProduct)
    npDotTimes.append(timeNpDot)
plt.plot(dotProductTimes, color = "blue", label = "dot_product()")
plt.plot(npDotTimes, color = "red", label = "np.dot()")
plt.title("Time spent computing dot products.")
plt.xlabel("vector size")
plt.ylabel("Time")
plt.legend()
plt.show()

Exercise 7. p-norms.

In mathematical analysis we often use the concept of normed spaces. Of particular interest are p-norms (also called $l_p$-norms) in $\mathbb{R}^n$, defined by

where $\textbf{x} = (x_0,...,x_{n-1})$ is a vector in $\mathbb{R}^n$. Particular examples of p-norms are the taxicab norm (also called Manhattan norm) with $p = 1$ and the usual Euclidian norm with $p = 2$. We can also define the maximum norm (also called the supremum norm or infinity norm) by

For each of these norms, we can of define unit circles by the equation \begin{equation} ||\textbf{x}||_p = 1. \end{equation}

Let us focus in the two-dimensional space $\mathbb{R}^2$. Plot in the same figure the unit circles corresponding to $p = 1, 2, 4, 16, \infty$.

Exercise 8. Shapes of arrays.

An important attribute of numpy arrays is their shape. Compute the shapes of the following arrays.

a = np.array([1,2,3,4])
b = np.array([[1,2],[3,4],[5,6]])
c = np.array([[[1,2,3,4], [5,6,7,8], [9,10,11,12]], [[13,14,15,16], [17,18,19,20], [21,22,23,24]]])

We now briefly discuss four nice features of numpy: vectorization, broadcasting, array indexing and boolean masks. These will prove invaluable for many numerical tasks we shall encounter later.

Vectorization

The numpy library was designed with support for array programming in mind. This means that numpy functions that apply to scalars can usually be used with arrays without the need of loops. Such functions are often said to be vectorized and usually result in performance gains when compared to loops. A typical example is the task of adding all the elements of an array. With a for loop, we would write something like

myArray = np.arange(1,101)
total = 0
for element in myArray:
    total += element

This can be compared with the simplicity of the np.sum() function:

total = np.sum(myArray)

As we have discussed in the beginning of this notebook, all familiar operations like +, *, /, //, ** are vectorized in numpy. But that is not all. Most mathematical functions in numpy are. The general process of vectorization works as

As an example, below we generate an array of 7 equally spaced points in a circumference and compute their cosines.

angles = np.linspace(0,2*np.pi, 7)
cosines = np.cos(angles)
print("angles: ", angles)
print("cosines: ", cosines)

Exercise 9. Vectorizing functions.

Often the functions we want to apply to arrays are functions we ourselves create. These are not vectorized by nature. The good news is that numpy has a commodity function np.vectorize() that vectorizes other functions. It works like this.

myVectorizedFunction = np.vectorize(myFunction)

Together with python's lambda expressions for defining anonymous functions, it can save us lines of code. As an example, the two pieces of code below are equivalent.

def myFunction(x):
    return x+3
myVectorizedFunction = np.vectorize(myFunction)
myArray = np.linspace(0,10,11)
result = myVectorizedFunction(myArray)
print("Result with standard function definition: ", result)

is equivalent to

myVectorizedFunction = np.vectorize(lambda x: x + 3)
myArray = np.linspace(0,10,11)
result = myVectorizedFunction(myArray)
print("Result with lambda expression: ", result)

Your task is to vectorize the bitwise NOT. We can think of this operation as the following map

Complete the missing steps below.

Broadcasting

A particularly useful feature of numpy arrays is broadcasting. This is the ability to operate with arrays of different dimensions. As a first example, how would we add the same number to all entries of an array?

myArray = np.array([1,2,3])
myInt = 10
myNewArray = myArray + myInt
print(myNewArray)

Numpy therefore understands our intention here: by adding a number to an array, we mean adding this number to all entries in the array. We thus optimize our code by avoiding an iteration. Let us consider a second example.

array1 = np.array([[1,2,3],[4,5,6]])
array2 = np.array([10,20,30])
mySum = array1 + array2
print("array1: ", array1)
print("array2: ", array2)
print("mySum = array1 + array2 ", mySum)

Again, we are operating on arrays of different shapes (array1 is (2,1), while array2 is (1,)). This is no problem to numpy, since it understands that the first dimension in array1 is equal to the first (and single) dimension in array2: both are 3. Now let us see an example of two arrays that cannot be broadcast together. If you run the code below, you should read an error message.

array1 = np.array([[1,2,3],[4,5,6]])
array3 = np.array([10,20])
mySum = array1 + array3

We see that numpy considers array3 as a line vector with two components. Since array1 can be seen as a 2X3 matrix, each line in array1 has three components. So we see what the problem is: we cannot add vectors with two components to vectors with three components. Pretty intuitive. But what if I wanted to add array3 to array1 columnwise? Could we just transpose array3 and achieve that? Let us test it.

mySum = array1 + np.transpose(array3)

No deal again! In order two operate on two arrays, it is required that they have compatible dimensions. From the documentation, two dimensions are compatible when

1. they are equal, or
2. one of them is 1.

We can understand what this means by checking the shape attribute of arrays array1, array2 and array3. Let's do this.

print("array1: ", array1, "Shape of array1: ", array1.shape)
print("array2: ", array2, "Shape of array2: ", array2.shape)
print("array3: ", array3, "Shape of array3: ", array3.shape)

Exercise 10. Add array3 to array1 columnwise.

How would you add array3 to array1 columnwise in the previous example? (Hint: you might want to look at np.reshape)

Exercise 11. Create the following arrays in numpy and guess their shape.

You can always check your answer using the shape attribute. \begin{equation} A = \begin{bmatrix} \begin{bmatrix} 1 & 2\\ 3 & 4 \end{bmatrix} & \begin{bmatrix} 5 & 6\\ 7 & 8 \end{bmatrix} \end{bmatrix} \end{equation}

Exercise 12. Broadcasting the arrays from the previous exercise.

If you created the arrays from the previous exercise correctly, it should be possible to broadcast $A$ to the shape of $B$ and $C$ to the shape of $D$. So numpy should be able to compute $B - A$ and $D - C$ for example. Perform these operations. Does the result agree with our intuition?

Array indexing

We often want to acces elements inside an array, be it for referencing or assignment. Numpy gives us several ways to achieve this. Before we actually get to array indexing, let us quickly review indexing in numpy. The canonical way to access an element of an array by index is just passing the element's indices separated by commas between square brackets. So for example, if we want to access the number 3 in the matrix $A$ below to change it to 9,

we just type

A = np.array([[1,2], [3,4]])
print("A before change: ", A)
A[1,0] = 9
print("A after change: ", A)

Exercise 13. Find the number 5.

Find the number 5 in each of the arrays in exercise 11. Change it to 42.

Numpy also supports slicing, a feature we already saw when studying lists. This allows us to capture many elements of an array without having to write their indices explicitly. Consider the array below.

Suppose we wanted to capture just a block of A, namely

We can easily do this with slicing.

A = np.arange(1,17).reshape(4,4)
B = A[0:2,0:2]
print('A: \n', A)
print('B: \n', B)

Note that, just as with lists, the lower limit of a slice is inclusive, while the upper limit is exclusive. None of this is new, since lists also have these capabilities. Let us introduce then array indexing, a feature of numpy that lists don't have. We consider the following example. We are given a one dimensional array of size ten and we want to extract the elements at positions 3 and 7. We can do the following.

myArray = np.arange(1,11)
indexArray = np.array([3,7])
newArray = myArray[indexArray]
print("myArray: ", myArray)
print("newArray: ", newArray)

Here we have used an array (indexArray) as an index to another (myArray). The name "array indexing" is thus quite literal. Indexing arrays with other arrays can be quite useful as the exercise below shows.

Exercise 14. Logarithmically spaced samples.

You are given a large one dimensional array whose values increase linearly with the array index. You want to take a sample from this array that is (approximatelly) logarithmicaly spaced. Code the function below that fulfills this task.

Exercise 15. Using array indexing with multidimensional arrays.

In this exercise we will experiment with indexing multidimensional arrays. The code snippet below generates a 10x10 array.

myArray = np.arange(0,100).reshape(10,10)
print("This is the shape of myArray: ", myArray.shape)
print("This is myArray:\n")
print(myArray)

Your goal is to use array indexing to generate the following subArray:

subArray = np.array([[myArray[3,3], myArray[3,4], myArray[3,5]], [myArray[4,3], myArray[4,4], myArray[4,5]], \
                    [myArray[5,3], myArray[5,4], myArray[5,5]]])
print("\n")
print("This is subArray: \n")
print(subArray)

Boolean masks

Numpy has a very powerful indexing tool based on boolean arrays called boolean masks. It is similar to vectorization, but now it is a boolean operator that is being vectorized. Assume we have a boolean function $b(x)$ that evaluates to True or False. The vectorization process then yields

Thus, we get as a result a boolean array with entries equal to true wherever $b$ evaluates to True in the given array, and False otherwise. As a simple example, run the code below.

myArray = np.array([1,2,3,4,5])
myBooleanArray = myArray > 3
print("myArray: ", myArray)
print("myBooleanArray: ", myBooleanArray)

The beautiful thing about boolean arrays is that they can be used for indexing. If we pass a boolean array as an index to a given array of same shape, the result is a new array consisting only of the elements where the boolean array is True. This is what we call a boolean mask. Using our previous example, let us create a new array from myArray applying the boolean mask myBooleanArray.

newArray = myArray[myBooleanArray]
print("newArray: ", newArray)

Exercise 16. Find the even elements of an array.

Write a function that takes an integer array as input and returns a boolean array that tells us where its even elements are.

Exercise 17. Logical operations with arrays.

The process of boolean masking can be seen as creating filters that we put over arrays so that only entries satisfying the filters' conditions pass through. Since these are actually logical filters, we can think of composing them using logical operations. Intuitively, we expect that something like this should work.

a = np.array([[True, False], [False, True]])
b = np.array([[True, True], [False, False]])
print("This is a: ", a)
print("This is b: ", b)

Then we could try

c = a and b
print("This is c: \n", c)

As you should see by yourself, this does not work. How would you achieve the same purpose? What about the logical or operator?

Other particularly useful logical operations built into numpy are the all, any and where functions. The meaning of these words in english should give you a good idea of what they do. In any case, an example will come in handy.

a = np.array([True, False, True])
b = np.array([False, False, False])
c = np.array([True, True, True])

print("This is a: ", a)
print("This is b: ", b)
print("This is c: ", c)

We now test where are the True values of a, whether any value of b is True, and whether all values of c are True.

print("Where are the True values of a? ", np.where(a))
print("Is there any True value in b? ", np.any(b))
print("Are all values in c True? ", np.all(c))

Exercise 18. Find common divisors.

You are given a very large array of integer numbers (call it intArray) and a tuple of positive integers (call it divisors). Your task is to find where in intArray are the numbers that are divided by all of those in divisors. Write a function that performs this task. Your function should return a dictionary whose keys are the indices of the desired elements of intArray and whose values should be these elements.

Exercise 19. Tic-tac-toe.

In this exercise we will write a program that plays tic-tac-toe with itself randomly. This will be a good opportunity for you to use many of the numpy array features we have seen up to this point. We set up the game as follows. The game starts with a 3x3 grid filled with zeros. The value 0 means that the corresponding spot in the grid is free. The game proceeds with player 1 placing a 1 in a free spot of her taste. Then player two places a 2 in a free spot of her own taste and so on. At every move, your program must check if the game must continue or not, either because there is no free spot anymore or one of the players has won. We have done some of the work for you. Fill in the missing steps.

The meshgrid.

In this section we are going to enhance our visualization skills with the help of the numpy meshgrid. Our goal here is to generate stuff like 3d plots, but you will find that the meshgrid is far more useful than that. These plots rely on specifying a function $z = f(x,y)$ at a region in the plane, so we have to be systematic on how we generate the points $(x,y)$ in the grid. The meshgrid does this for us.

Suppose you want to compute and visualize the values of a function $f(x,y)$ in the square $[0,4]$x$[0,4]$. You will have to generate a grid like this.

What the meshgrid does is generate this grid in a systematic way. It decomposes the matrix above as

Let's see how to do this in Python code.

xArray = np.arange(5)
yArray = np.arange(5)
print("This is xArray: ", xArray)
print("This is yArray: ", yArray)

Now we create our grid components.

xGrid, yGrid = np.meshgrid(xArray, yArray)
print("This is xGrid: \n", xGrid)
print("This is yGrid: \n", yGrid)

Now let us create the array in python corresponding to the paraboloid $f(x, y) = x^2 + y^2$.

zGrid = xGrid**2 + yGrid**2

Finally, this is how we plot the paraboloid in the square $[0,4]$x$[0,4]$.

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
newFig = plt.figure(figsize = (7,7))
ax = plt.axes(projection='3d')
ax.plot_surface(xGrid, yGrid, zGrid)
plt.title("This is my first 3d plot!")
ax.set_xlabel('This is the x axis')
ax.set_ylabel('This is the y axis')
ax.set_zlabel('This is the z axis')
plt.show()

Exercise 20. Practice with meshgrid.

Make a 3d plot of $f(x,y) = \sin(\pi^2 xy)$ in the square $[-1,1]$x$[-1,1]$.

Exercise 21. How related are two positive integers?

In this exercise we wish to visualize how related two positive integers are. We will use the following criterium. Given two positive integers $x$ and $y$, the degree to which they are related will be given by

where gcd stands for greatest common divisor and lcm stands for the least common multiple. Make a contour plot showing how th positive integers are related in the square $[1,100]$x$[1,100]$.